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Configurations of supercooled liquids residing in their local potential minimum (i.e. in their inherent structure, 
IS) were found to support a non-zero shear stress. This IS stress was attributed to the constraint to the energy 
minimization imposed by boundary conditions, which keep size and shape of the simulation cell fixed. In this 
paper we further investigate the influence of these boundary conditions on the IS stress. We investigate its 
importance for the computation of the low frequency shear modulus of a glass obtaining a consistent picture 
for the low- and high frequency shear moduli over the full temperature range. Hence, we find that the IS 
stress corresponds to a non-thermal contribution to the fluctuation term in the Born-Green expression. This 
leads to an unphysical divergence of the moduli in the low temperature limit if no proper correction for this 
term is applied. Furthermore, we clarify the IS stress dependence on the system size and put its origin on a 
more formal basis. 


I. INTRODUCTION 

Upon decreasing temperature supercooled liquids display 
a dramatic increase in their shear viscosity. This effect 
can be directly related to an increase of the relaxation 
time of shear stress fluctuations. In his seminal work 
Goldstein 1 pointed out the importance of the potential 
energy landscape over the N particle configuration for 
the slow relaxation dynamics of supercooled liquids. At 
a sufficiently low temperature the system is assumed to 
be close to a local minimum of this landscape (called an 
inherent structure). 

This picture naturally implies the presence of two time 
scales in supercooled liquids: a fast relaxation process 
associated with vibrational motion of the system around 
an inherent structure and a slow relaxation process corre¬ 
sponding to thermally activated hopping to another min¬ 
imum in the potential energy landscape accompanied by 
a rearrangement of a relatively small number of particles 
in the system. In a non-equilibrium situation applied 
strain might induce the disappearance of a local mini¬ 
mum facilitating such a rearrangement and leading to 
stress relaxation. While this fact has motivated investi¬ 
gations on the influence of strain on the local potential 
minim ePHU, several studies have implicitly made use of 
residual stresses in inherent structures present even at 
equilibrium in vest igating the magnitude 6 and the relax¬ 
ation dynamic^® of the inherent stresses as the system 
samples different minima. As pointed out by Abraham 
and Harrowell, the notion of a shear stress of an inherent 
structure is a priori far from obviousP if a system which 
can be brought arbitrarily close to a local energy mini¬ 
mum (say by an appropriate, numerical energy minimiza¬ 
tion of a simulation), one would naively expect that all 


global shear stresses in the system are zeroed by this min¬ 
imization procedure. However, it has been argued that 
boundary conditions impose a constraint on the energy 
minimization i.e. the shape and the size of the contain¬ 
ment of a supercooled liquid forbids to zero all stresses. 
In a computer simulation this constraint amounts to a 
particular choice for the simulation box. It is not clear 
in which sense this IS shear stress can be considered as a 
predecessor of the stresses supported by a deformed su¬ 
percooled liquid/glass or to which extent it determines 
the stress relaxation process in a non-equilibrium situa¬ 
tion at all. Therefore, a better understanding of the IS 
stress is highly desirable. A first discussion of it has also 
been given in [9[. Among other things, the authors found 
that the magnitude of the IS stress is surprisingly, essen¬ 
tially independent of temperature, scales with a certain 
power of the system density and with the inverse system 
size. 

The aim of this note is to lift the inherent structure 
stress on a more formal footing which will help us to 
further clarify its origin and to discuss in what sense it 
has an influence on computations of viscoelastic proper¬ 
ties. To this end, we will present a consistent picture for 
the proper calculation of low and high frequency shear 
moduli of glass-forming liquids. This work is organized 
as follows: In section II we begin our discussion by ap¬ 
plying the Irving-Kirkwood formula for IS configurations 
and tracing back the remaining stresses to the choice of 
boundary conditions in a formal sense. In section III, we 
proceed by identifying external mechanisms which bias 
the computation of shear moduli upon decreasing tem¬ 
perature in a supercooled liquid and provide calculations 
of these moduli for glass forming liquids. In section IV 
we summarize our results and conclude with a discussion 
of the physical meaning of the IS stresses. 
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II. THE IRVING-KIRKWOOD FORMULA FOR 
INHERENT STRUCTURE CONFIGURATIONS 

According to Irving and Kirkwood, the instantaneous 
stress tensor in a configuration with particle mass nii , 
positions restricted to a volume V is given 

1 n i N N 

a ctP = y 'y ' m i v i,a v i,P ~ ^ 1 ^ ' r ij,aFij,p i (1) 

2—1 2=1 _7 = 1 

where r.y = r ? ; — rj, Vi = r,; and fy is the pair forces 
exerted on particle i by particle j. The Greek indices 
refer to the cartesian component of the corresponding 
vector/tensor. In an IS configuration (i.e. in a mechan¬ 
ically stable packing) the particles are at rest and the 
inherent structure stress a IS is solely determined by the 
second part of ([l ). We assume periodic boundary condi¬ 
tions in a cubic box with side length L i.e. we redefine 
r ij,a = t i,a ~ fj,a + fiij,aL. Here, riij denotes the vector 
which minimizes the distance between particles i and j 
where its components are ±1 or 0. Inserting the peri¬ 
odic boundary conditions in the configurational part of 
([T]) leads to 
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By using Newton’s third law Fij = —Fji and renaming 
of indices, this can be rewritten as follows. 


is 

a aP 


N N N 


2V 




( 3 ) 


i=l j=l 


where Fj = Y^^=i Fij the net force acting on parti¬ 
cle i. This quantity vanishes by definition for an IS. In 
a simulation it is arbitrarily small in the sense that the 
magnitude of the net force on a particle is bounded by 
the smallest force tolerance at which the used minimiza¬ 
tion algorithm (e.g. a conjugate gradient solver) still con¬ 
verges. Therefore, the IS stress is approximately given by 

L N N 

a ct/3 ~ ~2V 53 ni 3 ,aFij,p ■ (4) 

*=1 j=l 

As we will discuss subsequently, expression Q has a 
straightforward physical interpretation. Note, that the 
cu-component of the vector n l:j is only nonzero if a parti¬ 
cle close to the a = L /2 boundary of the simulation box 
interacts with the periodic image of a particle residing 
at the opposite (a = —L/2) boundary or vice versa. For 
instance, if we choose a = x, /3 = y in a cartesian co¬ 
ordinate system, expression Q is nothing else than the 
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Figure 1. Variance of inherent structure stress in 2D (top, 
soft sphere system) and 3D (bottom, binary LJ system) both 
at temperature T=0.5 as a function of system size. Black 
circles indicate the full Irving-Kirkwood expression and red 
triangles the approximate boundary formula Q. The insets 
show the same data on a log-log scale with a linear fit which 
has a slope of —1. 


total force component in y-direction exerted by the right 
boundary layer of the simulation box on its left counter¬ 
part. This shows analytically in which sense the choice 
of boundary conditions (i.e. the shape of the simulation 
box) determines the IS stress. We test approximation 
© numerically by comparing the mean squared IS stress 
calculated from the full Irving-Kirkwood expression to 
the one calculated from approximation Q by performing 
molecular dynamics simulation of glass forming systems 
for different system sizes, temperatures and dimensions 
(see appendix A for details on the simulations). The to¬ 
tal inherent structure stress is very well approximated by 
equation © (see Fig.l). As it can be read off from Fig. 
1, the fluctuations of show a \/V system size depen¬ 
dence, which has already been described described in [ 5 ]. 
A heuristical explanation is given in appendix C based on 
equation ([4]). The investigation of the IS in amorphous 
materials dates back to ideas of Stillinger and Weber 11 
and has proven to have various applications including the 
investigation of rate processes in low-temperature amor¬ 
phous substances, the formulation of equations of state 
for supercooled liquids, macroscopic transport proper¬ 
ties etc. (see e.g. [T2].jT5] and El)- While equation Q 
seems to identify the IS stress to be a mere (negligible) 
boundary effect, we will discuss that it is felt through¬ 
out the system and has a major influence on macroscopic 
quantities. Our discussion will focus on the influence of 
IS stresses on the numerical computation of elastic con¬ 
stants, more specifically on the shear moduli of a glass 
forming material. 
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III. ELASTIC CONSTANTS 

The shear modulus, G, describes the response of a ma¬ 
terial to shear stress and is defined as the ratio between 
shear stress to shear strain. If a material is subjected 
to oscillatory shear deformation, the shear modulus is a 
function of the excitation frequency oj, i.e. G = G(w). 
Its high- and low frequency limits are characteristic for 
a material’s mechanical behavior: the infinite frequency 
shear modulus, Goo, describes the response to an instan¬ 
taneous, affine deformation. It is not directly measurable 
experimentally, since a deformation at truly infinite fre¬ 
quency cannot be applied in practice. Goo should not 
be confused with the experimentally reported high fre¬ 
quency modulus which always refers to the shear modulus 
at the highest obtainable frequencies.^ The temperature 
dependence of Goo is relatively weak and its value mainly 
depends on the microscopic details (atomistic potential) 
of the system. The zero frequency limit, Go, describes 
the ability of relaxing stresses on a long time scale. Since 
a liquid in equilibrium does not support any stresses, Go 
is zero for a liquid in equilibrium but finite for a solid. 
Therefore, the zero frequency modulus can be regarded 
as an indicator for solidity. It depends strongly on the 
thermodynamic state and is therefore sensitive to tem¬ 
perature changes as a material approaches its melting 
point. The situation is more complicated for glassform¬ 
ing materials as they do not show a sharp solidification 
transition. In the following, we will summarize and ex¬ 
tend previous results for Goo and Go of a supercooled 
liquid and discuss their behavior over the full tempera¬ 
ture range. Furthermore, we will investigate the contri¬ 
bution stemming from the IS and discuss in what sense 
it contributes to properties of the low temperature glass. 
Throughout this section the simulation results are re¬ 
ported for the two-dimensional soft sphere system with 
N = 512 particles (see appendix). 

A. General remarks 

The infinite frequency shear modulus is analytically given 
by the so-called Born-Green expression^ which is well 
defined and yields non-vanishing results in both the solid 
and the fluid phase. This expression is given by 


Goo — pk B TP 



for an isotropic system with the hydrostatic pressure P 
and a pair potential </>(r). Further utilizing isotropy and 
performing an orientational average this can be simplified 

tcjnnu 


for a two dimensional system. For a soft sphere system 
with a purely repulsive pair potential of the type r~ n 
equation <[6j) leads to 19 

TL — 2 

Goo = pk B T H-— (P - pk B T) . (7) 

The low frequency shear modulus is given by equat ion © 
corrected by the so-called fluctuation tern:P^o|2i], i.e. 



In the following we will provide an extensive discussion 
of these quantities for a glass forming liquid in different 
temperature regimes. 


B. High temperature liquid 

For temperatures well above the melting point of the 
system, the low frequency shear modulus vanishes, i.e. 
the system does not sustain non-zero stresses on a long 
time scale. In this situation the fluctuation term cancels 
the Green-Born expression and the high frequency shear 
modulus can be calculated by the shear stress fluctua¬ 
tions: 



Note that the ensemble average of the shear stress 
( a xy ) in equilibrium always vanishes. 

C. Supercooled phase 

Upon further cooling the system below its melting point, 
the supercooled regime is entered. The fluid is not in 
its true thermodynamic equilibrium being the crystalline 
phase but is said to be in a metastable equilibrium in the 
sense that time translational invariance holds and two 
time correlation functions (such as the stress-stress au¬ 
tocorrelation function, G(t) = (cT xy (t)a xy ( 0))) decay to 
zero within experimentally available time windows (see 
fig. 2). As the particle motion becomes increasingly slug¬ 
gish, relaxation becomes slower and happens on two time 
scales: vibrational degrees of freedom lead to a fast re¬ 
distribution of stresses (i.e. an initial decay of C(t) to 
the lowest possible value possible for a particular con¬ 
figuration in place). This process is followed by a slow 
relaxation associated with particle rearrangements on a 
mesoscopic scale. As the system is still (quasi-)ergodic in 
the sense that it finds a way to redistribute stresses such 
that G(f) fully decays to zero, the formulas for the elastic 
moduli (equations ([ 7 ]) and ©) are still valid. However, 
since the relaxation of stress correlations takes an increas¬ 
ingly long time scale, the liquid assumes a viscoelastic 
behavior and the low frequency shear modulus departs 
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from its zero value. According to equation ([8 ), this is as¬ 
sociated with a decrease in the shear stress fluctuations. 
It also means that equation © is not appropriate to com¬ 
pute the high frequency modulus anymore, but equation 
(5) has to be used. The fact that the common expression 
(9) looses its validity in the so lid phase has been pointed 
out by several authors L 9 I 19 I 22 I 


D. Glassy phase 


If the system is further cooled below the glass transition 
temperature T g a full relaxation of two time correlation 
functions cannot be observed anymore within the exper¬ 
imentally available time window. The calculation of the 
shear modulus of a material in this glassy state was ex¬ 
tensively discussed by Williams ! 22 * 23 ^ We briefly recap the 
physical picture considered by the authors there: The 
phase space of the system is divided into Njj subsys¬ 
tems. Every subsystem is in equilibrium but between 
the domains the system is out of equilibrium. The prob¬ 
ability density of the domain a is given by / a (T) = 
Sa(r) exp( ~| g(r)) with Z a = / dTs a (T) exp(—/37?(T)), 
where T is the phase space coordinate, H the Hamilto¬ 
nian of the system and s a a switching function which 
is equal to unity if T lies in the domain a and zero 
otherwise. The probability distribution of the entire 
system is given by a composition of the single-domain 
distribution weighted with a nonequilibrium weight, i.e. 

/( r ) = T,a=i w a.fa(T) and J2a=\ i w a = 1- Equation Q 

holds for the (equilibrium) subdomains only. The authors 
of this study further derived that the infinite frequency 
shear modulus of the system is given by 

V Nd 

Gp = GqqJ — ({<J xv ) n - (<Jxy) n ') , (10) 


where the subscript / denotes the rule to average the 
Green-Born expression over the distribution / and the 
subscript a means an equilibrium average over the do¬ 
main a. This makes an accurate and meaningful calcu¬ 
lation of the low frequency modulus for a glass sample a 
very subtle task as it would require to consider the single 
phase space domains. Under the assumptions that every 
simulation is sampling its own, single domain and that 
the set of prepared samples representatively reflects the 
distribution of the weights w a , one could estimate (10) by 
simple time averages. This approach was seemingly taken 
in j9j and in [TS]. However, as already pointed out in [23 
this method is very sensitive to the used time over which 
averages are taken. Therefore, we mention an alterna¬ 
tive approach, also presented in [22] : the low frequency 
modulus is given by the Green-Born expression minus the 
stress-stress autocorrelation function C{t) at t = 0 (equa¬ 
tion Q), which drops to a non-zero plateau value for a 
broad class of glass forming liquids. This means that the 
stress fluctuations relative to the frozen-in stresses in this 


domain are considered. The autocorrelation is now cal¬ 
culated up to a cutoff time t c which is much larger then 
the relaxation time of the fast processes in the sample. 
Finally, the fluctuation term in ([8])) is corrected by C(t c ): 

Go = G °°-^T (<**»> - W 2 - C '(*c)) • (U) 

It should be noted that this procedure has the advantage 
that it is not sensitive to the chosen cutoff since C(i) is 
almost constant for a broad time interval (meaning that 
ageing of the glass is negligible on the time scale of in¬ 
terest for the investigated model system). The physical 
picture behind this correction is the following: each in¬ 
dividual glass sample contains frozen-in stresses, which 
are induced due to the initial conditions and preparation 
procedure of the sample. For instance, a fast cooling 
protocol pushes a liquid out of equilibrium very rapidly. 
This does not leave enough time for stress relaxation pro¬ 
cesses to occur leading to significant stresses in the glass 
sample which might not be present if a slower cooling 
rate would have been used. While these residual stresses 
are sometimes deliberately introduced during the manu¬ 
facturing procesiP - and influence mechanical/rheological 
measurements on individual glassy materials, they are 
not a characteristic property of the material but a re¬ 
mainder of its production process. Correcting for the 
frozen-in stresses by subtracting the plateau value of the 
stress-stress correlation function removes this contribu¬ 
tion from the low frequency modulus such that Go re¬ 
mains a quantity which is characteristic for the material 
irrespectively of its history. As we will see in the next 
section, the IS stress may also bias the computation of 
the elastic moduli. Since, we want to investigate this 
effect in more detail, we focus on systems where correc¬ 
tions for the frozen-in stresses in the form of (|8j) play 
a minor role only. This is the case for the considered 
system. While the self-intermediate scattering function 
does not decay to zero anymore, when the liquid passes 
the glass transition temperature 2 -^, the stress-stress au¬ 
tocorrelation function C(t) decays to values which are 
negligible for a potential correction according to equa¬ 
tion ([8]) (see fig. 2). The physical reason for this behav¬ 
ior might be, that stress relaxation events happen only in 
local rare events, but the global stress-stress correlation 
decays since very few rearrangements on the boundaries 
of the simulation box lead to a large change in the IS 
stress contribution as discussed in section I. 


E. Low temperature limit 

The low temperature limit of the amorphous solid is of 
wide interest in the research community and subject to 
extensive scientific effort (for an overview see the cor¬ 
responding sections in 26)]. In the following, we will 
extend our discussion on the calculations of the elastic 
properties to the low temperature regime. We will use 
the present considerations about the shear moduli as a 
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C(t) 



expression has to be evaluated at the inherent structure 
configuration and the number in the subscript refers to 
the n-th derivative. Transforming to re-scaled coordi¬ 
nates r i}Cl — r™ a = (/cbT) 1 / 2 ?-' a and inserting (13) and 
(14) into equation © yields 


J S a fay + VkBT^iy, 1 ^ + •■■) 

exp -& 2 Sr ' r '^j ex P ^ b — &™r'r'r' — ...^ dr'/ 

J s Q exp ^ — ^$™ r 'r'^ exp (^—■^—^—$1 s r'r'r' — dr'. 

(15) 


Figure 2. Stress-stress autocorrelation function C(t) = 
Wxy(t)<Txy(0)) for high (T = 0.6, black) and low tempera¬ 
ture (T = 0.3, red). While at T = 0.6 the systems starts to 
develop a two step relaxation process, a pronounced plateau 
can be observed at T = 0.3. Note, that for the system under 
consideration C(t) drops to values close to zero even at tem¬ 
peratures below the nominated glass transition temperature 
T ~ 0.35. See text for an explanation. Numerical data have 
been fit by a stretched exponential aexp(— bt c ) + d for long 
time scales. 


tool to identify different mechanism contributing to the 
solidification process of a glassy material and clarifying 
the role of the IS stress played in it. 

For a subdomain a which is assumed to be in equilib¬ 
rium, we follow the calculation of Lutsko 27 and compute 
the shear stress fluctuation in the canonical ensemble in 
the low temperature limit. 

(°ly) a = Y a l drSa < e ^ (-j£r) ’ (12) 

where Z a = J drs a ex p and $ is the potential 

energy of the system. We expand both potential energy 
and stresses around its inherent structure, i.e., 


$ = *\t=T* 

1 < 9 2 $ 

2 dr-i^ a drj 


<9$ 


dr. 


in, a - G>)+ 


(n,a ~ r i%)( r i,P - TV j,) + ... , (13) 


and 


1 d 2 c 


r is da xy 
dri a 


2 dri^ a dr 


{n, a ~ G>)+ 

.is 

(n,a ~ ~ r{ J) + . 


(14) 


where we have used the summation convention for in¬ 
dices. In the following, we will use 4> 7,s and a™ n as 
a short notation, where the superscript means that the 


Expanding the second exponential factor in powers of 
k B T (in both denominator and enumerator), all integrals 
are Gaussian integrals, which are solvable analytically. 
Note that odd moments of these Gaussian integrals van¬ 
ish. Finally, expanding the quotient in equation (15) in 
powers of k B T leads directly to the following result for 
the stress-stress fluctuation in the low temperature limit: 


V 


k B T 




xy/f 


to ?y a ? 

k B T 


Wn 


+ A + Bk B T + 0((k B T) 2 ) 


(16) 


where a™: a is the inherent structure contribution from 

x y 

particles in subdomain a and the subscript / has the same 
meaning as in equation (( |10[ )). At this point we have to 
make further assumptions in order to estimate the low 
temperature limit of the fluctuation term. Again we are 
left with the problem of identifying the different subdo¬ 
mains in order to estimate the weights w a and to sum 
over (a 2 y ) a accordingly. We will make the previously 
mentioned assumption that the prepared samples reflect 
the distribution of these weights and that each simula¬ 
tion predominantly samples its own single domain such 
that the total inherent structure stress of one simulation, 
a™ approximates a™ ,a . Therefore, the low temperature 
limit of the fluctuation term is estimated by ensemble 
averaging over statistically independant starting config¬ 
urations. Hence, the quantity A is the linear term in the 
perturbation expansion and given by 


/ dcr™ da™ \ 

y Iql j/3 J 


((rLr'e)) 


1 IS ( d<T iy \ 

3 xy \dr ia drjpdr kl dr u J 


(( r'iar'jpr'^r'ig )) 


+ a 


is 

xy 


( ^j S y \ 

y dr ln dr jB J 


{(r' a r' j0 ))- (17) 


The contributions to the linear term B are discussed 
later. We have introduced the Gauss bracket of a 
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function g which is defined as follows ((g(r'))) = 

1 J dr’g(r') exp (— $ 2 Sr ' r ') > where c is a normalization 
constant such that ((!)) = 1 holds. Note, that 


{{ r iafjp)) — ((*^2 ) 


-1 


(18) 


where the right hand side is nothing but the inverse of the 
Hessian matrix. The Hessian matrix has zero eigenvalues 
due to the fact the the net force on the system is zero 
and is therfore not invertible. A commonly used solution 
to this issue is, physically holding one particle fixed (i.e. 
excluding it from the sums in all calculations), which has 
no effect on the free energy.^ All higher moments like 
{{ri a Tj 0 r kl ru)) can be traced back to (18) using Wick’s 
theorem, e.g.: 

(( r ia rjpr kl ri$)) = (( r ia rjp))((ri n ris)) + 

((f’ial'ky')') ((^j finis')') “k ((y'ial'is')') ((Tj 7)) ■ (19) 


We find that the first term in equation ( |16| ) is exactly 
the inherent structure contribution to the stress fluctu¬ 
ations. Since these fluctuations are essentially tempera¬ 
ture independent, this term would lead to a divergence 
of the low frequency shear modulus. This means that, 
under the constraint of the given containment (or the 
boundary conditions of the simulation box) the particles 
possibly cannot be packed in a way which zeros the total 
stress as discussed in the previous section. We note that 
this constraint is imposed on the entire box so it does 
not affect the results of local elastic quantities where the 
subvolume is embedded in a larger system (e.g. see [29] or 
EDI). At high temperatures this term also does not play 
an important role: for formal reasons due to the 1/T 
prefactor, for physical reasons due to permanent stress 
redistributions activated by thermal motion occurring at 
high temperatures. These stress fluctuations are usually 
much larger than its IS contribution (see figure ([3])). We 
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identify the IS contribution to the shear stress fluctua¬ 
tion as geometric frustration of the system, accompany¬ 
ing the kinematic frustration, which has been discussed in 
the previous section on the nonequilibrium glassy state. 
Again, this contribution is not an inherent material prop¬ 
erty but reflects external constraints on the relaxation 
dynamics becoming significant at low temperatures. The 
second term in ( |16[ ) is temperature independent and the 
benefit of equation (17| is that, we can read off how the 
stress fluctuations obtained from a simulation should be 
corrected to extract the true zero temperature properties 
of the material when frozen-in stresses can be neglected. 
Alternatively, the material properties can be calculated 
from the first term in which does not contain the 
IS shear stress itself but only its derivative. The latter 
is nonzero irrespectively of any constraint on the system. 
This term has already been deduced in [27], however it 
has so far not been numerically tested for amorphous sys¬ 
tems to the best of our knowledge. Even though equa¬ 
tion (16) is strictly speaking only valid in equilibrium 


at low temperatures, we propose to correct the zero fre¬ 
quency shear modulus even in the glassy regime, since 
the inherent structure stress plays a predominant role al¬ 
ready at this temperature range as can be seen in figure 
[3] The second term, B, in equation ([T6| describes the 
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slope at which the fluctuation term departs from its zero 
temperature limit. Among all terms contributing to this 
order of (fcsT) we neglect those which contain IS con¬ 
tributions (see Appendix B) and obtain the temperature 
dependence of the fluctuation term close to T = 0. 

As a conclusion, we have obtained a detailed picture of 
the computation of the high and low frequency shear 
moduli of a glass forming system: While in the high tem¬ 
perature regime the fluctuation term cancels the infinite 
frequency shear modulus, it decreases in the supercooled 
regime. This means the the glass forming material de¬ 
velops an elastic behavior. Due to the non-thermal IS 
contribution the fluctuation term would increase again 
upon further cooling in the glassy phase. As we correct 
for this contribution, we observe a small decrease of the 
fluctuation term towards the low temperature limit (see 
figure [4]). The difference between G M and the fluctua¬ 
tion term is depicted in figure 5. At high temperatures 
Go is essentially zero meaning that the fluid does not 
support stresses on a long time scale. Entering the su¬ 
percooled regime the material behaves elastically, which 
is mirrored by an increase of the zero frequency modu¬ 
lus. As the temperatures is reduced further this elastic 
behavior becomes more pronounced (see figure [5]). 


IV. CONCLUSION 


Figure 3. Fraction of stress variance stemming from the 
IS contribution, increasing over 50% when the system enters 
the supercooled regime. At this point stress redistributions 
through thermal activation become so slow that the exter¬ 
nal constraint, set by the boundary conditions becomes non- 
negligible and requires a correction to the shear moduli. 


In this note we discussed the origin of the inherent struc¬ 
ture stress and traced it back to the particular choice 
of boundary conditions. In view of equation @ the IS 
stress can be understood as the net forces penetrating 
the boundary of the simulation box. Notably, the for- 
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Figure 4. High frequency shear modulus (blue, open circles) 
according to equation 0- Filled, black circles show simula¬ 
tion results for the fluctuation term from ensemble averaging 
with the proper IS corrections applied (see text). The red 
square is the zero temperature limit according to the first 
term in equation <©• The red line corresponds to perturba¬ 
tion theory results for the fluctuation term of the Born-Green 
theory. The low frequency shear modulus is given by the dif¬ 
ference between black and blue data points. We distinguish 
four temperature regimes: Regime I is the high temperature 
regime at which equations 0 and 0 hold equally to calcu¬ 
late the high frequency modulus. In regime II the system is 
supercooled. The fluctuation does not fully cancel Goo any¬ 
more meaning that the material develops an elastic behavior. 
Regime III is the glassy state at which both frozen-in and 
IS stresses bias the calculation of the fluctuation term and 
are corrected to extract inherent material properties irrespec¬ 
tively of external constraints or history dependence. Regime 
IV is the low temperature glass phase at which the fluctuation 
term can be computed with equation (161 



Figure 5. Zero frequency shear modulus obtained from sub¬ 
tracting the fluctuation term from the high frequency shear 
modulus according to equation |8|. Black circles show simu¬ 
lation results with proper IS corrections applied and the red 
square is the zero temperature limit according to the first 
term in equation (161. The red line corresponds to perturba¬ 
tion theory results for the fluctuation term of the Born-Green 
theory. 


mula Q exactly coincides with the so-called “Method of 
Planes” definition of the s tress evaluated at the bound¬ 
aries of the simulation box. 31 This analytic expression 
allowed us to explain the scaling of the IS stress with 
the system size (see appendix). Moreover, it might serve 
as starting point for understanding other properties of 
the IS stress, e.g. its temperature independence or its 
absolute magnitude. Finally, we would like to address 
the question which physical meaning lies behind the IS 
stress. It is obvious that in a sufficiently low tempera¬ 
ture regime the configurational part of the stress tensor 
becomes dominant over the kinetic contribution. Start¬ 
ing from the Green-Kubo formula for the viscosit}®, 
= FiT Jq 50 (t)cr a ,y (0))dt , it has been noted that the 
IS stress autocorrelation describes the onset of highly vis¬ 
cous behavior as a liquid enters its supercooled regimeP 
Hence, as the system explores the potential energy land¬ 
scape, hopping from one minimum to another, the un¬ 
derlying IS stress fluctuations determine its viscosity in 
the linear response regime. In view of equation Q this 
means that a rearrangement of forces at the boundary 
of the system contains enough information to character¬ 
ize the inherent structure which the system currently re¬ 
sides in. However, two caveats seem to come along with 
this notion. First, the IS stress itself does not provide 
any obvious information on which time scale its auto¬ 
correlation will decay or at which frequency the hopping 
between the energy minima occurs. It is merely associ¬ 
ated with a geometric frustration which is imposed on 
the system by choosing particular boundary conditions. 
Secondly, one should not conclude from Q that cooper¬ 
ative rearrangements taking place at the transition from 
one minimum to another primarily happen at the bound¬ 
ary of the system. They occur anywhere in the system 
but the energy minimization zeroes all net forces on par¬ 
ticles locally, leading to a redistribution of the boundary 
forces. In this sense the boundary conditions introduce 
a frustration which is present everywhere in the system. 
In view of the previous discussion the following, over¬ 
all picture for the role of the IS stress emerges: the IS 
stresses and the redistributions of forces might amount to 
a macroscopic, scalar observable which characterizes the 
inherent structure such that its autocorrelation mirrors 
the path of the system through its configuration space, 
but does not contribute to the viscosity via any associ¬ 
ated timescale as it does not contain dynamical informa¬ 
tion itself. However, it affects the elastic properties of 
the system as it leads to a diverging contribution to the 
fluctuation term of the Born-Green theory. The reason 
for this behavior is that the IS contribution to the fluc¬ 
tuation term is of a non-thermal origin. If this term is 
divided by fcgT, it leads to a divergence of the fluctua¬ 
tion term and therefore to results for the elastic moduli 
which are biased by the external constraint imposed on 
the system. Since every rheological (or mechanical) mea¬ 
surement goes together with a change of the boundary 
of the system and as there is no experimental indication 
for a low-temperature divergence of the infinite frequency 











shear moduli in glassy systemiP^LUJ ^ see ms appropriate 
to remove the IS contribution to the shear modulus in or¬ 
der to correct for external constraints in the preparation 
procedure of the IS. 
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Appendix A: Model systems and simulation details 


The chosen model system and simulation details es¬ 
sentially correspond to those used by Abraham and 
HarrowelP investigating glass forming model systems in 
2 TW and 3lM Molecular dynamics simulations were 
carried out at fixed volume, temperature and parti¬ 
cle number, using the LAMMPS packaged We used 
equimolar, binary mixtures of particles interacting via a 
purely repulsive potential in 2D, 0y(r) = e (—) , with 

e = 1, (Tn = 1, (T 12 = 1-2,(722 = 1-4 and a Lennard- 
Jones potential in 3D, <f>ij{r) = 4 e ((^-) 12 — (^ i ) -6 ) 
with e = 1, (Jn = 1 , (j 12 = 1.1,(722 = 1.2. A cut¬ 
off r c was used at 4.5crn and the potential were shifted 
such that they vanished at r c . The mass of the particles 
were mi = 1 and m 2 = 2 and the used densities were 
p = 0.747 in 2D and p = 0.75 in 3D. Both systems were 
simulated at T = 0.5 which is below the freezing tem¬ 
perature of both systems and the particle numbers were 
chosen to be IV = 512,1024,2048,4096,8192 in 2D and 
N = 1024,2048,4096,8192 in 3D. Reduced units are used 
for length L* = L/an, temperature T* = fcgT/e, energy, 
E* = E/e and pressure P* = Pa 3 /e. Throughout all 
simulations, periodic boundary conditions are applied. 


Appendix B: Perturbation calculations 


In the following we list the ter ms w hich contribute to 
the linear order term in equation (16 1 , B = JA =1 -Si- 


Si — n Nii/,ai cr x!/,a2 l ^a3a4(i5«6(( r ai r Q2)) \ \ J_ ]_ r a, 
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IF 
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IF 


B'J Q&xy,cti&xy,Oi20L3Ci4 


if 


Note that the calculation of this first order term is a 
somewhat lengthy and complicated task: Firstly, higher 
derivatives of the total potential $ and of the Irving- 
Kirkwood expression of the stress are required. However, 
these can be computed in a straightforward way for sim¬ 
ple pair potentials (like the purely repulsive soft sphere 
potential used in this study) by using software which is 
capable of performing symbolic computations. Secondly, 
some of these objects (e.g. the fourth derivative of the 
potential ^ ai a 2 a 3 a 4 ) contains a very large number of en¬ 
tries. In the present study we investigate systems with 
512 particles in two dimensions meaning that this tensor 
has approximately (10 3 ) 4 = 10 12 entries. Summing over 
this number of entries and even storing the object itself 
seems to be an intractable task. However, the majority 
of entries are zero, since the derivative with respect to 
three (or more) different particle coordinates vanishes for 
a pair potential. Additionally, the sequence of the per¬ 
formed derivatives does not influence its value according 
to Schwartz’s theorem. A third problem arises due the 
use of Wick’s theorem: the term B$ requires the calcula¬ 
tion of the eighth moment. Splitting this in all possible 
pair combinations yields in total 105 terms but the prob¬ 
lem can be simplified using symmetry properties of the 
involved objects. In the following we provide further de¬ 
tails on how one of the terms is calculated in practice. 
The other terms can be handled in a completely anal¬ 
ogous way. For instance the term H 4 can be split as 
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follows: 

B4 @xy,a 2 ^asa^asae 

{{AA A) {AA A A A) + 

{AAA) {AAAAA) + A ■ ( B1 ) 


peaked around a characteristic distance, it can be esti¬ 
mated that ro ~ a/p 1 ^ in d dimensions. Therefore, we 
assume the interaction between particle i and j to be 
determined solely by density of the system but to be in¬ 
dependent of its size. Summing over the index j in Q, 
we obtain 


The first term of the sum cancels B±. Contract¬ 
ing over ai and introducing the abbreviation T as = 
0 xy,ai (( r ai r a 3 ))i this can be rewritten 


rs 


L 

V 


N b 

EC ' 1 

i=1 


(Cl) 


Ba = — — a 
24 ‘ 


Xy , Oi 2 ^0:30:40:50=6 


( t «3 {AAAAA) + 


T a4 «r' 


OL 5 016/ 

Aa 3 A r A) + -) ■ (B2) 


Since the sequence of the derivatives of $ does not matter 
all terms in the brackets are the same which leads to 


-®4 — 24 ^a 3 a4a 3 ae (^a 3 {{ r a 2 r a i r a 5 r a e ))) ‘ 

(B3) 

Employing Wick’s theorem another time, i.e.: 

{AAAAA) = {AAa , e )) {AAA) + 
{AAA) (AAA) + AAA) {AAA) ( B4 ) 

and using the permutability of indices again, we conclude: 


£4 = ~T a3 T ai <f> a3aia5a6 {AAA) • ( B5 ) 

This expression can be handled numerically bearing in 
mind that only those entries of < & a3a4a6Q , 6 are non-zero for 
which at least three particle indices are identical. With 
these simplifications the first order contributions B\ to 
E ?7 can be calculated. We note that, even though the 
slope of temperature dependence is in good agreement 
with the medium temperature data the sample to sample 
fluctuations of the first order term B is much larger than 
those of the zeroth order term A. 


Appendix C: Volume dependence 


where F x y L is the y-component of the net force act¬ 
ing on a particle i for which Xi < L/ 2 exerted by the 
periodic images of the particles residing in the vicinity 
of the Xj = —L/ 2 boundary. Nb is the number of par¬ 
ticles in the boundary layer of the simulation box i.e. 
( Xi,yi ) € [L/2 — r c ,L/2] x [—L/2,L/2], where r c is the 
interaction range of the potential. We also multiplied by 
a factor of two to account for the reversed situation where 
particle i resides at Xj > — L/2. Clearly, the average over 
different configurations vanishes, (cr^y) ~ 0 but we find 
for the average magnitude of the IS stress: 


<«) 2 >»4(£A 


N h 


V 2 


N b \ 

? x=L /2 ST' t?x=L /2 

,y 2—/ * j,v 

3 =1 / 


(C2) 


The sum in (C21 extends over particles which are in close 


vicinity of each other. Therefore, the forces F X ~ L ^ 2 can¬ 
not be regarded to be uncorrelated but might be assumed 
to be weakly anti-correlated over distance of a few par¬ 
ticle diameters. From a physical point of view, this is 
perfectly sensible, as will be made clear by the following 
intuitive argument. Assume, a particle i\ is subjected 
to a strong, positive force = a > 0). Then, a 

particle i 2 , which is not to far away, should push with 
the same force Ff ~A 2 ~ —a in the reverse direction in 
order to maintain the mechanical stability in which the 
inherent structure configuration resides per definition. It 
is clearly not possible that all forces perfectly cancel in 
that manner. We denote the number of unbalanced forces 
by N c which should scale as the number of particles in the 
corner of the simulation box i.e. N c ~ r 2 L d ~ 2 . Applying 


this argument to equation ( C21 leads to 


In this appendix we provide an intuitive explanation 
for the scaling behavior of the IS stress variance. We start 
by briefly recapitulating a simple density scaling argu¬ 
ment for plastic flow of amorphous solids made in m and 
[55] , which was used previously to explain the density 
scaling of the IS stress. We consider a purely repulsive 
pair potential <p(r) = ( a/r) n and the force magnitude be¬ 
tween particle i and j is given by Fij = d< A' J ^ ~ r A~ l ' 
If we denote the probability distribution of distances 
r at density p by p(r, p) , the mean distance is given 
r-o(p) = f rp(r,p)dr. Assuming that p(r, p) is strongly 


<(<) 2 > 


V 2 


2 / N b N b \ 

E f x=l/2 e p x = L / 2 


\ i=1 


3 =1 
N c 


V 2 


E( Kv L/2 ) ) * ( C3 ) 


w=i 


where the remaining forces in (C3) are uncorrelated. 
Invoking the central limit theorem to the presumably in¬ 
dependent remaining forces results in the following seal- 






10 


ing estimate 



(C4) 


which interestingly reveals the IS stress fluctuations 
to be a V~ l effect although shown to arise from a 
contribution originating from a constraint imposed on 
the boundary of the box. This scaling behavior coincides 
with the one numerically found by the authors of [5]. 
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